# 利用专家知识
"""
在特征工程中，专家知识非常重要。虽然机器学习旨在减少对专家规则的依赖，但这并不意味着要忽视领域知识。
领域专家可以帮助识别出比原始数据更有信息量的特征。例如，在预测机票价格时，原始数据包括价格、日期、航空公司、出发地和目的地等信息。
机器学习模型虽然可以从这些数据中构建出不错的模型，但可能无法捕捉到某些关键因素，比如度假高峰期和假日期间的机票价格通常会更高。
虽然一些假日的日期是固定的（如圣诞节），但其他假日的日期可能取决于月相（如光明节和复活节）或官方规定（如学校放假）。
如果只使用公历记录日期，模型可能无法学习到这些事件的影响。通过添加一个特征来编码航班是否在公休假日或学校假期之前、之中或之后，
可以将先验知识编码到特征中，辅助机器学习算法。即使最终发现假日信息对机票价格没有显著影响，添加这一特征也不会有害处。
"""
import mglearn
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.preprocessing import OneHotEncoder, PolynomialFeatures

'''
关于预测 Andreas 家门口的自行车出租数量的任务，这是一个典型的机器学习问题，可以通过以下步骤来解决：
数据收集：
 使用 Citi Bike 提供的公开数据集，这些数据通常包括每次租车的时间、地点、持续时间等信息。
 确保数据集包含 Andreas 家门口的自行车租赁站的数据。
特征工程：
 时间特征：提取日期、时间、星期几、是否为节假日等特征。
 天气特征：添加天气数据，如温度、降雨量、风速等，因为这些因素会影响自行车的使用。
 历史数据：考虑过去几天或几周的租车数据，作为特征。
 专家知识：利用常识，比如工作日和周末的租车需求可能不同，节假日的租车需求可能更高。
模型选择：
 选择合适的机器学习模型，如线性回归、随机森林、梯度提升机等。
 使用交叉验证来评估模型的性能，确保模型具有良好的泛化能力。
模型训练和评估：
 使用训练集训练模型，并在验证集上评估模型的性能。
 调整模型参数，优化模型性能。
预测和部署：
 使用训练好的模型对给定日期和时间的租车需求进行预测。
 将模型部署到实际应用中，实时预测 Andreas 家门口的自行车出租数量。
通过这些步骤，可以有效地预测 Andreas 家门口的自行车出租数量，帮助他提前了解是否有自行车可供使用。
'''

# 我们首先将这个站点 2015 年 8 月的数据加载为一个 pandas 数据框。我们将数据重新采样
# 为每 3 小时一个数据，以得到每一天的主要趋势：

citibike = mglearn.datasets.load_citibike()
print("Citi Bike data:\n{}".format(citibike.head()))


# 下面这个示例给出了整个月租车数量的可视化
plt.figure(figsize=(10, 3))
xticks = pd.date_range(start=citibike.index.min(), end=citibike.index.max(),
 freq='D')
plt.xticks(xticks, xticks.strftime("%a %m-%d"), rotation=90, ha="left")
plt.plot(citibike, linewidth=1)
plt.xlabel("Date")
plt.ylabel("Rentals")

'''
观察数据可以发现，每24小时的白天和夜间有明显区别，工作日和周末的模式也不同。对于时间序列预测任务，通常用过去的数据预测未来，划分训练集和测试集时，
以某个特定日期为界，之前的数据作为训练集，之后的作为测试集。这里用前184个数据点（23天）作为训练集，剩余64个数据点（8天）作为测试集。
预测任务的输入特征是日期和时间，输出是接下来3小时内的租车数量。计算机存储日期常用POSIX时间，即从1970年1月1日00:00:00起至现在的总秒数。
可以尝试用这个单一整数特征作为数据表示。
'''

# 提取目标值（租车数量）
y = citibike.values
print(citibike.head())
# 利用"%s"将时间转换为POSIX时间
X = citibike.index.astype("int").reshape(-1, 1)

# 我们首先定义一个函数，它可以将数据划分为训练集和测试集，构建模型并将结果可视化

# 使用前184个数据点用于训练，剩余的数据点用于测试
n_train = 184
# 对给定特征集上的回归进行评估和作图的函数
def eval_on_features(features, target, regressor):
 # 将给定特征划分为训练集和测试集
 X_train, X_test = features[:n_train], features[n_train:]
 # 同样划分目标数组
 y_train, y_test = target[:n_train], target[n_train:]
 regressor.fit(X_train, y_train)
 print("Test-set R^2: {:.2f}".format(regressor.score(X_test, y_test)))
 y_pred = regressor.predict(X_test)
 y_pred_train = regressor.predict(X_train)
 plt.figure(figsize=(10, 3))
 plt.xticks(range(0, len(X), 8), xticks.strftime("%a %m-%d"), rotation=90,
            ha="left")
 plt.plot(range(n_train), y_train, label="train")
 plt.plot(range(n_train, len(y_test) + n_train), y_test, '-', label="test")
 plt.plot(range(n_train), y_pred_train, '--', label="prediction train")
 plt.plot(range(n_train, len(y_test) + n_train), y_pred, '--',
          label="prediction test")
 plt.legend(loc=(1.01, 0))
 plt.xlabel("Date")
 plt.ylabel("Rentals")

# 我们之前看到，随机森林需要很少的数据预处理，因此它似乎很适合作为第一个模型。我
# 们使用 POSIX 时间特征 X，并将随机森林回归传入我们的 eval_on_features 函数。
from sklearn.ensemble import RandomForestRegressor
regressor = RandomForestRegressor(n_estimators=100, random_state=0)
plt.figure()
eval_on_features(X, y, regressor)

'''
在训练集上的预测结果相当好，这符合随机森林通常的表现。但对于测试集来说，预测结
果是一条常数直线。R2 为 -0.04，说明我们什么都没有学到。发生了什么？
问题在于特征和随机森林的组合。测试集中 POSIX 时间特征的值超出了训练集中特征取值
的范围：测试集中数据点的时间戳要晚于训练集中的所有数据点。树以及随机森林无法外
推（extrapolate）到训练集之外的特征范围。结果就是模型只能预测训练集中最近数据点
的目标值，即最后一次观测到数据的时间。
'''

'''
显然，我们可以做得更好。这就是我们的“专家知识”的用武之地。通过观察训练数据中
的租车数量图像，我们发现两个因素似乎非常重要：一天内的时间与一周的星期几。因此
我们来添加这两个特征。我们从 POSIX 时间中学不到任何东西，所以删掉这个特征。首
先，我们仅使用每天的时刻。如下所示，现在的预测结果对一周内的每天都具有相同
的模式：
'''
X_hour = citibike.index.hour.reshape(-1,1)
eval_on_features(X_hour,y,regressor)

# R2 已经好多了，但预测结果显然没有抓住每周的模式。下面我们还添加一周的星期几作为特征

X_hour_week = np.hstack([citibike.index.dayofweek.reshape(-1, 1),
 citibike.index.hour.reshape(-1, 1)])
eval_on_features(X_hour_week, y, regressor)

'''
现在我们的模型通过考虑一周的星期几和一天内的时间捕捉到了周期性的行为。它的 R2
为 0.84，预测性能相当好。模型学到的内容可能是 8 月前 23 天中星期几与时刻每种组合
的平均租车数量。这实际上不需要像随机森林这样复杂的模型，所以我们尝试一个更简单
的模型——LinearRegression
'''

eval_on_features(X_hour_week, y, LinearRegression())

'''
LinearRegression 的效果差得多，而且周期性模式看起来很奇怪。其原因在于我们用整数
编码一周的星期几和一天内的时间，它们被解释为连续变量。因此，线性模型只能学到关
于每天时间的线性函数——它学到的是，时间越晚，租车数量越多。但实际模式比这要复
杂得多。我们可以通过将整数解释为分类变量（用 OneHotEncoder 进行变换）来获取这种
模式
'''

enc = OneHotEncoder()
X_hour_week_onehot = enc.fit_transform(X_hour_week).toarray()
eval_on_features(X_hour_week_onehot, y, Ridge())

'''
它给出了比连续特征编码好得多的匹配。现在线性模型为一周内的每天都学到了一个系
数，为一天内的每个时刻都学到了一个系数。也就是说，一周七天共享“一天内每个时
刻”的模式。
利用交互特征，我们可以让模型为星期几和时刻的每一种组合学到一个系数
'''
poly_transformer = PolynomialFeatures(degree=2, interaction_only=True,
 include_bias=False)
X_hour_week_onehot_poly = poly_transformer.fit_transform(X_hour_week_onehot)
lr = Ridge()
eval_on_features(X_hour_week_onehot_poly, y, lr)

'''
这一变换最终得到一个性能与随机森林类似的模型。这个模型的一大优点是，可以很清楚
地看到学到的内容：对每个星期几和时刻的交互项学到了一个系数。我们可以将模型学到
的系数作图，而这对于随机森林来说是不可能的。
'''

# 首先，为时刻和星期几特征创建特征名称：
hour = ["%02d:00" % i for i in range(0, 24, 3)]
day = ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"]
features = day + hour

'''
然后，利用 get_feature_names 方法对 PolynomialFeatures 提取的所有交互特征进行命名，
并仅保留系数不为零的那些特征：
'''

features_poly = poly_transformer.get_feature_names(features)
features_nonzero = np.array(features_poly)[lr.coef_ != 0]
coef_nonzero = lr.coef_[lr.coef_ != 0]

# 下面将线性模型学到的系数可视化

plt.figure(figsize=(15, 2))
plt.plot(coef_nonzero, 'o')
plt.xticks(np.arange(len(coef_nonzero)), features_nonzero, rotation=90)
plt.xlabel("Feature name")
plt.ylabel("Feature magnitude")

plt.show()